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1. QUARTERLY REVIEW AND ABSTRACT 


FOREWORD 


Project Whirlwind 

Project Whirlwind at the Massachusetts Institute of Technology Digital 
Computer Laboratory is sponsored by the Office of Naval Research under Contract 
NSonbO. The objectives of the Project are (l) the application of an electronic 
digital computer of large capacity and very high speed (Whirlwind I) to problems 
in mathematics, science, engineering, simulation, and control, and (2) the study 
and development of component reliability in Whirlwind I. 


The Whirlwind I Computer 

Whirlwind 1 is of the high-speed electronic digital type, in which quantities 
are represented as discrete numbers, and complex problems are solved by the 
repeated use of fundamental arithmetic and logical (i.e., control or selection) 
operations Computations art executed by fractional-mirrosecond pulses in elec¬ 
tronic circuits, of which the principal ones are (l)the flip-flop, a circuit containing 
two v,»i uurn tubes so connected that one tube or the other is conducting, but not 
both. (2) the gate or coincidence circuit. ( 3) the magnetic-core memory, in which 
binary digits are stored as one of two directions of magnetic flux within ferro¬ 
magnetic cores. 

Whirlwind I uses numbers of lb binary digits (equivalent to about 5 decimal 
digits). I hi $ length was selected to limit the machine to a practical size, but it 
permits the computation of many simulation problems. Calculations requiring 
greater number length are handled by the use of multiple-length numbers. Rapid- 
access magnetic-core memory has a capacity of 32,768 binary digits. Present 
speed of the computer is 40, 000 single-address operations per second, equivalent 
to about 20,000 multiplications per second. This speed is higher than general 
scientific computation demands at the present state of the art, but is needed for 
control and simulation studies. 


During the past quarter S2 problem* made u*e >1 the computer time allotted to the 
Scientific and Engineering Computations (SfcLC) Group The computer logging desc ribed 
previously is now being utilised to produ. • .» biweekly report d time used by each problem, 
the processing is done entirely by WWI New routines lor displaying numerical results on 
the oscilloscope were introduced into the comprehensive system of service routines. 

Greater emphasis inow being placed on reliability and maintenance of WWI to 
provide the best possible performance. Increased reliability resulted in 95% usable appli¬ 
cations time comparril with an average of 9 3% Since Septcmber 1963. To this end a biweekly 
analyst . is being made of the types of interruptions occurring most frequently and those 
resulting in greatest loss of time. 

There were 32 students representing 12 Ml'I groups cm oiled in the DCL program¬ 
ming course given this quarter Topics included were relative addresses, temporary 
storage, floating addresses, preset parameters, programmed arithmetic, cycle counters, 
buffer storage, automatic output, post morterns, and multipass conversion. The course 
on machine computation. Machine-Aided Analysis, had an enrollment of 69, primarily 
seniors in Electrical Engine. ring. Practice problems were done by the students on a KEAC 
and on WWI (simulating the three-address compute r developed for the 1964 summer session). 
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2. MATHEMATICS, CODING, AND APPLICATIONS 


2.1 Introduction 

During the period covered by tlu» report Si problems made use of the computer time 
allotted to the Scientific and Engineering Computation (SliEC) Group. Progress reports as 
submitted by the various programmers arc presented in numerical order in Section £. Z. 
Of these, 21 (189. 191, £03, £10-£18, ZZO-ZZZ, 225, and 228*232) represent new problems 
that are being described for the first time. Fourteen problems (107, 159, 166, 184, 190, 
£00, £08, £10, £11, £1£, £11, £15, ££0, and £££) have been completed The results of £4 
problems will be included in academic theses (one S. B. , two S. M. , 18 PI.. D. or Sc. D). The 
results of 11 problems have been or will be submitted for publication in technical journals. 

The computer logging described in Summary Report No. 37 (page 1£) is now being 
utilized to produce a biweekly report of time used by each problem with the processing be¬ 
ing done entirely by WWl. 

New routines for displaying numerical results on the oscilloscope (for photographic 
output) have been introduced into the comprehensive system of service routines. A descrip¬ 
tion of the routines will also be found under Problem 100. 

2.2 Problems Being Solved 

100 COMPREHENSIVE SYSTEM OF SERVICE ROUTINES 

Th. comprehensive ,y,tem of service routine, has been developed by the Scientific 
and Engineering Computation (SlEC) Group to aimplify the process of coding for WWI. The 
system now in use. called CS II, was described in Summary Reports No. 36. 37. 38. and 39. 

Since the reader will find references in some of the reports below to the number 
system used in CS U. the following brief description is included here for the reader's con¬ 
venience. 

(m.n) numbers shall mean numbers which are of the form z = x • 2 y where x is an 
m-binary-digtt number and y is an n-binary-digit number. For example. (24.6) signifies 
a two-register floating-point system dealing with numbers of 24 significant binary digits 
(roughly seven decimal digits) with magnitudes between 2 63 and 2‘ 64 . 

Arithmetic involving these (m,n) numbers is carried out by means of (m.n) inter- 
pretivc subroutines . These subroutine, enable the programmer to write coded pro^T^I 
using (m.n) number, a. easily as. or even more easily than, he might write program, in 
the single-length fixed-point (14,0) number system which is bulli into Whiriwind l. 


Autom atic Biweekly 

During each computer session a logging tape is automatically punched out under the 
control of the computer. This logging tape contains the following information; (1) the time 
one problem began and ended, (£) whether the operator erased before rcad-ln, (3) the num¬ 
ber of the tape read in (this number also identifies the problem, programmer, and kind of 
tape), (4) whether a tape was converted, (5) whether a postmortem was given, (6) time used 
for testing terminal equipment, and (7) the computer's "down" time. The time is either in 
houes and minutes or standard time. 

A routine has been written which will read in these logging tapes and (1) sum up all 
the time used by each individual problem, (£) compute how many problems were run, (3) 
compute how many programs were run, (4) compute the time used for drum testing, .mag¬ 
netic-tape testing, calibration of the scopes, photoelectric-tape-reader check, etc., (5) 
compute the lost time, (6) sum up the over-ail time used, and (7) compute the percentage 
of usable computer time. 

These results arc combined with special information, stored on the auxiliary drum, 
about the problem titles and prescribed biweekly report formal to produce a "delayed" 
punched paper tape. This paper tape is then used to produce on Multilith the SfcEC Group's 
biweekly report. 

Even though the present system is working satislactorily, further development will 
continue to make the system faster and more flexible. 

Automatic Oscilloscope Output 

Automatic-output instructions for displaying numerical results on the scope have 
been added to the CS 11 vocabulary. The instructions SOA and iSOA (which stand for Scope 
Output Alplianumcrics, the 1 indicating that an interpretive mode has been selected) are used 
exactly like TOA and iTOA. Thus iTOA + l£3.45s means "type the contents of the MRAon 
the Flexowritcr as a signed number having three digits to the left of the decimal point and 
two to the right and follow it by a space" while iSOA + l£3. 45s will display the same number 
on the scope followed by a space. The "machine functions," carriage return, tab, and 
space, have much the same meaning on the scope as on the Flexowriter. A frame has room 
for 36 lines of numbers, and there arc four tab positions and room for 63 characters across 
the face of the scope. 

Some instructions are peculiar to the scope. The FRAME (iFRAME) instruction ad¬ 
vances the film in the camera, and, if SOA (iSOA) is used, it resets the display coordinates 
to the upper left-hand corner of the scope. The carriage return also does this if the number 
of lines exceeds 36. The COLUMN (iCOLUMN) instruction moves the "left hand margin 
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( MIT Project on Machine Methods of Computation) 


Table £-U. Current Problems Arranged According to the Mathematics Involved 
(* MIT Project on Machine Methods of Computation) 
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stop" to the next tab position to the right and resets the display position to the top line. If 
the next tab position is off the right-hand side of the scope, it acts like FRAME. 

106 mit seismic project 

A 9 discussed in previous Summary Report;.. Problem 106 is concerned with the 
invc stigatiun ol tiie use uf stall stical-analysis techniques in seismic - record interpretations, 
and in particular in the separation ol "rellections" from background interference on these 
records. More complete desc riptions of tile problem and the approaches used arecontained 
in "Detection of Reflections on Seismic Records by Linear Operators" (Wadsworth, Robin¬ 
son, Bryan, and Hurlcy--Geophysics, Vol. 18, No. 1, July 1953). This work is being 
carried out by the MIT Geophysical Analysis Group under the direction of Dr. S. M. Simp¬ 
son, 

In Summary Report No. 39 a list was presented of the various categories of com¬ 
putation in which Problem 106 has been interested. During the past 3 months emphasis 
has been put on three of these categories 

I. Computations designed to test techniques for probing the statistical and fre¬ 
quency behavior of a given input seismogram; 

3. Computations aiding the development ol the relationships between disc rete tinea r 
operators and linear electric filters; 

3. Computations designed to simulate and study various statistical noise models. 

Under 1 we have programmed a finely spaced, multiterm correlation cosine trans¬ 
form program, using a tallormadc programmed arithmetic to solve the speed-accuracy 
problem, and we have completed computation yielding a table for use in estimating signal- 
to-noisc ratios as functions of frequency from mulUlrace signal noise mixtures. 

Under l we have programmed the computation of inverse linear operators, the 
computation of linear operators from the location of their poles in the "z" plane, and the 
computation of linear operators from their real frequency power spectra. This last pro¬ 
gram, which is a means ol bypassing least-squares methods by spectral factorization, is 
not yet completed In addition, we have computed tables of linear operators with special 
characteristics: one set with high pass filter characteristics and another with narrow band 
pass characteristics. 

Under 3 we have computed a 1300-term stationary series by moving summation 
from a wavelet and a random-number tabic. From this scries, estimates of variance, 
correlation, and power spectra have been obtained for various sample lengths as an ex¬ 
perimental determination of variability. the results to be compared with observed seismic 
variability in these quantities. 
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The most important results obtained during this period have been, verification of a 
new concept of noise on actual seismograms; programs to exploit this concept, the tabulated 
ope.atois mentioned above; experimental information about seismic-wavelet stability; and 
experimental information about variability in estimates. 

For the future we plan to pursue these topics and to complete and use the spectrum- 
factorization program. 

I 07 (a) AUTOCORRELATION AND (b) FOURIER TRANSFORM 

INTEGRAL EVALUATION 

Routines developed under this problem by D. T. Ross of the MIT Servomechanisms 
Laboratory have been used by F. Raichlen of the MIT Hydrodynamics Laboratory ip the 
investigation of turbulent-velocity fluctuations in open-channel flow measured by means of 
a Pitot-tube pressure-cell combination. This project is sponsored by the Office of Naval 
Research under Contract No. NS-ori-07874. In the future, this project will be reported on 
under Problem III. 

During the past quarter it was found necessary to delete the secondary diaphragm 
of the transducer to increasr the mechanical stability of the system. This removal neces¬ 
sitated an increase in thickness of the primary diaphragm of the gage to bring the natural 
frequency of the system back to an acceptable value. A number of runs in a high-velocity 
flume were then made to investigate the effect of this change on the previous results ob¬ 
tained over the past year from the Digital Computei Laboratory. 

As before, the correlation curve obtained showed superimposed periodicities of 
approximately 180-cps frequency. These curves, however, were much cleaner in appear¬ 
ance as far as low-frequency components were concerned. This made possible a different 
approach to the computation of the macroscale of turbulence L x - By the application of a 
three-point smoothing process to the results obtained from the Digital Computer Laboratory 
the high-frequency periodicities were averaged out. The low-frequency component was 

removed from the curve by subtracting a cosine curve which had been fitted to this pre- 

—1( r 

dominant periodicity. A theoretical curve of R j - e was then fitted to the smoothed 
experimental correlation curve. From this the value of L x was obtained and found to be 
approximately the same as tnat computed previously by more laborious methods. 

108 AN INTERPRETIVE PROGRAM 

During the last quarter an effort has been started by Dr. J. H. Laning of the MIT 
Instrumentation laboratory to rewrite the algebraic interpretive program in its entirety. 
The object of this effort is to increase both the speed and versatility of the algebraic rou¬ 
tines. It is expected that an interim revised program will be available shortly in a form 
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suited to future development And extension*. 

Among the anticipated features of the new program are the following: 

1. The operation speed will be increased by an estimated factor of 5 to 10; 

2. The usable computer storage capacity will be considerably enlarged; 

3. The basic vocabulary of the program will be considerably increased: for ex¬ 
ample. function subroutines will now be designated by familiar words making the program¬ 
ming task somewhat easier. In addition, a number of new features will be included such 
as the automatic evaluation of sums, products, and integrals. 

At present the main logical structure of the program is about 75% complete. It is 
expected that this will be operating by the end of the present quarterly period. Following 
this, future revisions will be added as time and manpower permit. 

120 THERMODYNAMIC AND DYNAMIC EFFECTS OF WATER INJECTION INTO HIGH- 

TEMPERATURE. HIGH-VELOCITY GAS STREAMS 

This problem i. connected with the development of a potential gas-turbine compon¬ 
ent, called an "aerothermopressor." in which a net rise in stagnation pressure of a hot 
gas Stream is brought about by the evaporation of liquid water injected into a high-velocity 
region of the flow. The concepts underlying the operation of the aerothermopressor are an 
outgrowth of comparatively recent work in the field of gas dynamics, and its intended func¬ 
tion in the gas-turbine cycle is analogous to that of the condenser in a steam power plant. 

The device consists of a converging nozzle which accelerates the exhaust gases 
from the turbine into a circular duct of varying diameter terminated by a conventional 
conical diffuser, which recovers the kinetic energy of the flow before discharging it to the 
atmosphere. At the entrance of the duct (evaporation section), special Injectors deliver 
minute jets of water which are in turn atomized by the rapidly moving gas stream. The 
changes in state within the aerothermopressor are brought about by the simultaneous ther¬ 
modynamic and dynamic effects of (I) evaporation of the liquid water, (2) momentum and 
energy interactions between the phases, (j) friction, and (4) variations in cross-sectional 
area of the duct. Under proper circumstances, these effects bring about a net rise in 
stagnation pressure across the device. Further descriptions of this device may be found 
in earlier reports, beginning with Summary Report No. 32, Fourth Quarter 1952. 

The role of Whirlwind I in the successful development of the aerothermopressor is 
intimately connected with the determination o( performance characteristics of the device 
under all conditions of operation by means of a comprehensive one-dimensional analysis of 
the process. The principal objectives of this analytical study, which began late in 1952, 
are (I) to establish a firm theoretical foundation for an understanding of the process and 
(2) to obtain design data for experimental work. 
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During the past quarter, most of the effort was directed toward the planning and 
preparation of a completely revised Whirlwind program. This revision, which is currently 
about one-half completed, was motivated by obsolescence of the original program, a natural 
consequence of both expanded knowledge concerning the aerothermopressor and changes in 
the Whirlwind computer. The basic objective of the new program is to provide a rapid and 
convenient scheme of analysis, and the essential changes in the program include: 

1. Use of the fourth-order Runge-Kutta method for numerical solution of the dif¬ 
ferential equations, rather than the elementary Euler method; 

2. Change of independent variable from fraction evaporated (humidity) to axial 
distance along duct; 

3. Extensive use of the auxiliary magnetic drum for both computations andprogram 
organization, the latter in an attempt to avoid a plethora of punched tapes; 

4. Elimination of all interpretive instructions except those mandatory for arith¬ 
metical operations; 

5. Careful revision of all portions of the program involving iteration, for example, 
the subroutine used for determining roots; (The revised program is more than twice as 
fast as the original.) 

6 . Elimination of necessity for hand computation of a portion of the initial data by 
incorporation into program; 

7. Provision for the optional use of a concentric circular plug within the evapora¬ 
tion section of the aerothermopressor; 

8 . Revised treatment of the liquid temperature and liquid velocity to avoid instability 
toward the end of the calculation. 

Following completion of this program, it is planned to resume computations on a 
broad scale. Considerable experimental data from the medium-scale test facility at the 
MIT Gas Turbine Laboratory is being accumulated and awaits corroboration. Further, 
several procedures for calculating optimum performance of the aerothermopressor await 
exploitation. 

The aerothermopressor development program is being carried out at MIT under 
the sponsorship of the Office of Naval Research and is directed by Professor Ascher H. 
Shapiro of the Department of Mechanical Engineering. The theoretical aspects of the prob¬ 
lem treated by Whirlwind 1 are being carried out by Dr. Bruce D. Gavril. 

122 COULOMB WAVE FUNCTIONS 

A program hasbeen written by J. Uretsky and M. Rotenbergof the MIT Physics De¬ 
partment to generate Coulomb wave functions by the method of Brett, Yost, and Wheeler. 

I. Physical Review 49 . 175 (1936) 
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Runs are now being made to determine the region of usefulness for this method. 

The program has proven to be logically correct. However, there is some doubt ae 
to the convergence of the formulae being uaed. 

123 EARTH RESISTIVITY INTERPRETATION 

This problem is concerned with a means ol analysing earth-resistivity data. The 
procedure iato assume a fixed numberof homogeneous horizontal layers, together with the 
thicknesses and conductivities of these layers, and then to relax the Slichter kernel calcu¬ 
lated therefrom to the observed kernel by systematically changing the assumed parameters 
(see Summary Report No. 57). 

To investigate the behavior of the solution, several kinds of data have been run. 
Kernels were calculated theoretically for cases in which the second layer was much thinner 
than the firstorin which the conductivity contrast between successive layers was very small. 
These were rounded off to the accuracy of field data (larger random errors were added to 
some) before analysis. Kernels calculated for a number of layers greater than assumed in 
the analysis were also analyzed. In the statement of the problem, it is assumed that the 
conductivities of the first and last layers are known, as is the case. Data in which an incor¬ 
rect value of bottom-layer conductivity was used were analyzed. Also, an actual observed 
kernel, from an area in which drilling information was available for comparison, was run. 

At the same time the program was sped up so that an analysis takes about 4 minutes 
instead of 15 minutes, and several features were incorporated to eliminate the unstable 
oscillation of the solutions which sometimes results from an increase in the parameter- 
change multiplier. 

The behavior of the analyses in general seems satisfactory. However, if the initially 
assumed values of the parameters being sought are too far wrong, no solution will be ob¬ 
tained. because one of the denominators in the Crout calculation of the parameter changes 
is the result of subtraction of two nearly equal numbers. A change in the number system 
used from (24. 6) to (25. 5) or (26, 4) might help in some such cases. Provided that such is 
not the case, the analysis yields solutions which seem reasonable Intuitively and whose ker¬ 
nels match the initial kernels to a degree depending on the situation. Three-layer kernels 
fit to better than the accuracy of the data. Four-layer kernels, in the three-layer analysis, 
yield a "smoothed" solution. The field case fitted the drill information extremely well. A 
few more difficult cases and three more field cases are being run. 

The kernels for the field cases were integrated by Whirlwind from field-resistivity 
data taken in mining areas in the Southwest. 

This work is being carried out by K. Vozoff of the Department of Geology and Geo¬ 
physics under the sponsorship of the Project for Machine Methods of Computation and Nu¬ 


14 


RELEASE. CASE 06-1104 


MATHEMATICS. CODING. AND APPLICATIONS 


merical Analysis. 

126 A DATA-REDUCTION PROGRAM 

Problem 126 la a .ery large data-reduction program for uae In the Servomechanlams 
Laboratory. The over-all problem ia compoaed of many component aectiona which will be 
developed aeparately and then combined at a later date. Early efiorta were focused on the 
development of utility-type programs. These programs, which have been described in 
previous Quarterly Reports, include a fully automatic program to fit polynomials to arbi¬ 
trary empirical functions. a general-purpose Lagrange interpolation program, and a flexible 
and fairly elaborate post-mortem routine. The programs are being developed by Douglas 
T. Ross and William M. Wolf with the assistance of Miss Dorothy A. Hamilton. Servo¬ 
mechanisms Laboratory staff members. 

The basic data-reduction program has been completely checked using the mistake 
diagnosis routine (MDR.) to extract over 30 intermediate results for accuracy and logical 
checks. Work on this program has been dormant pending the development of the manual- 
intervention program (M1V). The M1V program is a general routine for the interpretation 
of manual interventions into the computer via insertion registers and activate buttons. It 
has facilities for several types of alpha-numeric and graphical scope outputs, sounding of 
audible alarms, and elaborate program modification. The MIV program will operate in 
conjunction with the basic data-reduction program (or any similar program), and the intent 
is to have the flexibility normally associated only with analog computers available lor 
instantaneous modification and monitoring of this class of digital-computer programs. 

The mistake diagnosis routine has been modified so that it now will operate from 
and on any drum group or combination of groups. The method of specifying parameters 
for break points has also been simplified and made more foolproof. A slight modification 
operates in conjunction with the MIV program to give scope outputs of MDR results as they 
occur, 

A "scope input" routine has been written which gives a type of two-dimensional 
analog input to the Whirlwind computer. The principle of the program is similar to the 
"flying spot scanners" used on analog computers. The equipment used is a 16-inch oscil¬ 
loscope under the control of the computer with a photocell mounted so that its field of view 
is the entire scope face. Then by programming a flying spot and asking whether or not the 
photocell "saw" the spot, the program can be made to follow an opaque pointer as it is 
moved in a random fashion over the face of the scope. Since the program displays the spot 
by digital coordinates, the tracking of the pointer constitutes analog input to the computer. 

Present work includes testing and use of the MIV program with the present basic 
data-reduction program and also extensive modification of the basic program to Include 
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several new equations. 

131 SPECIAL. PROBLEMS (STAFF TRAINING, DEMONSTRATIONS. ETC.) 

A payroll program ha*» been prepared for use m demonstrations. It computes hours 
worked, regular earnings, premium earnings, gross earnings, social security tax, with¬ 
holding tax, and net pay. 

There are two inputs to the program. One is a file of employee records containing 
employee name, numoer, rate of pay, and number of dependents. The other input, which 
the demonstration audiences will be invited to prepare, represents the weekly time cards. 

The direct printer is used to prepare pay roll checks and deduction slips on standard 

forms. 

During the past J months, eight groups visited the Laboratory. The affiliations of 
these groups are given in Section 5. 3. 

132 SUBROUTINES FOR THE NUMERICALLY CONTROLLED MILLING MACHINE 

Use of the computer during this quarter by J. H. Runyon of the MIT Servomechan¬ 
isms Laboratory has been restricted mainly tothe checking of milling-machine tapes. Two 
short subroutines for use with a feedrate-control routine were written and tested. A test 
program for this routine is being written. Programming for milling-machine tape prepar¬ 
ation for a two-dimensional cam has begun. 


memory in a form such that LSR MA b may be used to find the solution. The order of the 
polynomial, the number of points in the table of data, the location of the data on the drum, 
and the location of the matrix in core memory arc all program parameters Calculations 
are done in ( 30 - j • j) arithmetic. 

I 11 SELF-CONSISTENT MOLECULAR ORBITAL 

An improved version of the linear combination of atomic orbitals for self-consistent 
field method has been written by Dr. A. Meckicr •»( the MIT Solid State and Molecular The¬ 
ory Group. The pattern of the program is much like that indicated in Summary Report 
No. 3S but with three major modifications 

The new program does not presume the use of an orthonormalized set of orbitals as 
the basis of the calculation. The equations now have the following form for a closed-shell 
system: 

The molecular orbitals arc expressed as 


< M X >= l V x > v 


where v^(x) is a prescribed orbital and c^ j is a coefficient to be determined. In general, 


/v^(x)v x , (x)dx = A 


M' 


1 11 SA.EC SUBROUTINE STUDY 

A subroutine for solving simultaneous equations by Crout's method has been included 
in the library under the title LSR MA b, This subroutine will compute A~*B, where B is a 
matrix consisting of one or several columns. Check columns are provided to indicate the 
accuracy of the results. 

A new sine-cosine evaluation subroutine has also been included in the library. This 
routine, LSR FU 4b, differs from the older LSR FU 4 sin-cos routine in that the reduction 
of the given angle to the appropriate range is accomplished more rapidly, thereby reducing 
computing time in cases where many large angles arc involved. 

An improved version of LSR MA 5, Symmetric Matrix Inversion and Square-Root 
Inversion, ha9 been written, tested, and placed in the library. 

A subroutine, LSR MA 7, has been added to the library. Its function is to generate 
the augmented matrix of a set of simultaneous equations, the solution of which is the set of 
coefficients of the polynomial which is the best fit in the least-squares sense to some data 
given in tabular form. The data is stored initially on the drum. The matrix is left in core 


a nondiagonal but positive-definite matrix. The side condition that the (p's be orthonormal¬ 
ized is equivalent to 

C*A C = l . 

With p = CC* , the energy expression is 

E= Tr^Fi 1/iS Pfj>i G P(j\ ■ 

X a 

F is the one-electron energy matrix and 

G^ = i/t/“(l)^I)g u X<>(2)a(i)dx 1 dx 2 - I/2i/°(l)<r(l)g u X a (^)^i)dx 1 dx 2 - 1/2 Ju a (l) 

X (l)g u d <, U)|i( 2 ) d *| d * i 

with the restriction that the integrals be real. 

The conditions on p are now 

Tr/A = n = half the number of electrons 


IL 


17 







APPROVED FOR 


MATHEMATICS. CODING. AND APPLICATIONS - 

<^A )* * />A • 

The requirement that the energy be stationary with respect U infinitesimal variations of p 
leads to the necessary and sufficient conditions 

H/j A=Ap H 

wh ' re H X<r- F X<r *£, A** 0 

The substitutions H=A"* ^ HA* * ^ 

yield the equivalent equations 

H p = p H 

~L _ - 

p - p Tr p * n 

These are solved as were the set for the old program. That is, H is diagonalized, the rect¬ 
angular matrix C is formed from its n lowest eigenvectors, and /T = CC*. This is followed 
by a transformation to p which is used to form a new H, and so on in the cycle. Thus, ex¬ 
cept for a preliminary calculation ofA‘* 1 or.d intermediate transformations, the iteration 
pattern is like the one previously programmed. 

The second improvement in the program is the test for self-consistency. The first 
program tested the variation of total energy from itcration to iteration. This is numerically 
bad. The numbers representing the eigenvectors are stabilized long before that number 
which represents the total energy. In fact, there is enough roundoff error in a computation 
of the large total energy so that it is impossible to judge the important last decimal places. 
Self-consistency can be missed. 

In the new program, the considered quantity is 

Tr l p > -p) 2 

where p is the input matrix and p • is the output matrix for one cycle. 

Now 

Tr( p' -p )~ - Tr {p‘ -p )£( p<. p )£ 

STr(p' -p) 2 [Tr A] 2 

the inequality following from the positive.definiteness of A . The equations satisfied by 
any p imply 
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T r( p' • p ) 2 s TrS p - -TrpS^i + 8 pp ■ 

II 8p were really a first variation, the idempotency of p would insist that 

8p p 8 p * 8pp 

Every p has the same trace, so that 


Tr hp ■ 


0 . 


Therefore, if the initial and final matrices did differ by only a first variation, Tt( p‘ - p )* 
would be equal to -Tr8/> and would vanish. The smallness of this trace or, in fact, the 
smallness of Tr(^>' -p )*" is taken as a test for self-consistency. 

The third improvement in the program is an over-all one engendered by machine 
developments and programming maturity. The use of buffers, for example, should increase 
the accuracy; the working is more direct. Better use of storage has allowed bigger systems 
to be handled if a symmetry type is of order n., then 


n.fnj + 1) 


S 99. 


155 SYNOPTIC CLIMATOLOGY 

Our basic objectives as put forth in Summary Report No. 39 have not changed suf¬ 
ficiently to require a restatement of them here. 

Over the past quarter we h^ve made great strides in expanding the size matrix that 
can be solved. This in turn has permitted us to add to our initial set of parameters others 
which we feel contain a sufficient amount of information with regard to our problem. The 
present program for solving matrices allows us to go as high as an n x n matrix, where n 
is restricted only by the amount of storage available. Thus n is restricted to 145, using 
only the amount of storage ten auxiliary-drum groups afford. The matrix-solution technique 
which we make use of is that developed by Crout, and the program is only applicable to sym¬ 
metrical matrices. 

Up to the present time we have solved a matrix as high as 84 x 84. This was accom¬ 
plished in a little over 12 minutes, even though portions of the program had to bt. done in the 
interpretive mode. A check was made on the amount of roundoff error accumulated. It 
appears from this check that the results, with input data containing seven digits, are ac¬ 
curate to seven places. 

Our aim during the next quarter will be to Incorporate magnetic-tape storage, which 
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should make the solution of any size matrix possible. Our interest in attempting to solve 
higher and higher matrices is due to the fact that until now we have as yet to reach a point 
of diminishing returns, insofar as the number of variables used and the number of degrees 
of freedom lost are concerned. 

Programming is being performed by Mr. Herbert J. Brun and Mr. Robert G. Miller 
under the supervision of Professor Thomas F. Malone of the MIT Meteorology Department. 

159 water use in a hydroelectric system 

This problem, undertaken by J.D.C. Little of the MIT Physics Department, has 
been described in Summary Report No. 37. In the past quarter some final runs have been 
made. The work on this problem has been reported in an MIT Physics Department doctor¬ 
al thesis. A paper is now being prepared for submission to the Journal of the Operations 
Research Society of America. 

162 determination of phase shifts from experimental cross sections 

A phase-shift analysis is being made of the elastic scattering of protons by 0*^ by 
Dr. F. J. Eppling ol the MIT Laboratory for Nuclear Science. The phase shifts required 
to provide a fit to a number of experimentai angular distributions will be determined from 
a partial wave expansion of the differential cross section. 

During the past quarter, a program designed to find for any combinations of phase 
shifts the best fit to the experimental cross sections has been worked on. Several tapes 
have been run to test different parts of the program. It is expected that soon the entire pro¬ 
gram will be ready to be run for any energy desired. 

166 CONSTRUCTION AND TESTING OF A DELTA-WING FLUTTER MODEL 

The basic procedures employed in this problem for designing a low-aspect-rationing 
flutter model to simulate a given set of flexibility influence coefficients of a full - scale wing 
were described in Summary Reports 37, 38, and 39 by S. I. Gravitz and M M. Chen of the 
MIT Aero-Elastic and Structures Research Laboratory. 

The basic computations involved were to yield solutions for a set of simultaneous 
nonlinear algebraic equations for the problem considered. The computation of a three-bay 
static model simulating the flexibility influence coefficients of an actual delta-wing airplane 
has been completed, and the results show good agreement with experimental data. Such a 
flutter model is still under construction and will be tested in the near future. 

The computation was satisfactorily concluded early in this quarter. The basic com¬ 
putational procedures appeared in an S. M thesis submitted by S.I. Gravitz to the MIT Aero- 
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nautical Engineering Department. A more detailed report concerning the designing, con¬ 
struction, and testing, etc., will be published under Contract Noas53-712-C, Bureau of 
Aeronautics, Department of the Navy. 

167 TRANSIENT EFFECTS IN DISTILLATION 

Work ha* continued on the problems of unsteady-state behavior of distillation sys¬ 
tems which have been described in preceding Summary Reports. The major emphasis has 
been on utilization of results obtained previously using Whirlwind I to understand the behav¬ 
ior of the systems. Additional calculations were made when the analysis justified them, 
but no new programs were necessary. The bulk of the work was on batch distillation, and 
it has been quite fruitful. Programs used were those for equilibration of a system with 
holdup and for product takeoff assuming holdup finite or zero. 

First the batch-distillation operation was studied as a whole, considering the equil¬ 
ibration period as well as the takeoff period when holdup was present. Results were studied 
on the basis of the ratio of product to charge, at constant values of product composition. 
Calculations were made changing the holdup-to-charge ratio and using various combinations 
of "apor consumption \n equilibration and takeoff. The first conclusion was that to obtain 
any specified product at any holdup-to-charge ratio there was an optimum combination of 
vapor consumption in equilibration and takeoff which resulted in minimizing the total vapor 
consumption. Second, the minimum consumption was less, as the holdup was less. Con¬ 
sequently this conclusion would dictate the minimization of the holdup in a practical system 
(subject to the necessity of having some holdup). 

These results may also be used to predict approximately the optimum vapor con¬ 
sumptions for other operations of distillation systems with holdup. A graphical method for 
predicting the total vapor consumption has been developed. The optimum *'apor consump¬ 
tion in equilibration can also be estimated. In general, it appears that a column should be 
equilibrated at total reflux at least until the top-vapor composition equals the desired prod¬ 
uct composition. How much more vapor should be used in equilibration depends on the vari¬ 
ables of the system. However, only with very high vapor consumptions during takeoff was 
more than twice this amount of vapor required for equilibration. 

Another phase of the investigation was concerned with studying how the vapor con¬ 
sumption during equilibration affected the top-vapor composition as a function of these sys¬ 
tem variables: holdup-to-charge ratio, charge composition, relative volatility, and num¬ 
ber of theoretical plates. Using Whirlwind, data were obtained over a broad range of these 
variables. Using these data .long with analytical developments based on linear equilibrium 
relationships, a substantial understanding of this nonlinear system has been obtained. A 
method has been found for predicting the top-vapor composition approximately, as a func- 
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lion of vapor throughout. This ia based primarily on the initial rates of equilibration of 
various parts of the system. Consequently it can be calculated from the initial conditions. 
This method may well be applicable in understanding and predicting the transient behavior 
of continuous distillation systems. These are the major developments. They are being 
written up by J. F. O'Donnell of the MIT Chemical Engineering Department in his Sc.D. 
thesis. 


172 OVERLAP INTEGRALS OF MOLECULAR AND CRYSTAL PHYSICS 


Since the last report it has been necessary to rewrite and extend several of the pro¬ 
grams used in the calculations. In particular, production-style programs were written and 
tested for solving two types of secular equations involving real symmetric matrices. These 
equations are the ordinary variety (Type I), 

n n 


Z H. .Xjm = 

j=i l J 


j=i 


rv E. , 

jm mm' 


and the more general variety (Type II), 

n 


z 

j=l 


H,.X 
ij jm 


n 

I S..X. E 
j s j jm mm' 


where n is an integer from 1 to 32 and H and S are the given real symmetric matrices (S 

positive definite). The results are displayed photographically and consist of the E , the 
. mm 

eigenvalues, and the X. , the eigenvectors, subject to the normalizations. 


Type I. 


I", X i/ X Jk = 8 ik and Type II, 


z 

j=l 


x fi s ij x jk 


’Ik- 


These programs are available lor the use ol any interested persons. 

A second program which was written and tested is a subroutine for punching out pro¬ 
gram-specified sections of the auxiliary magnetic drum as an edited, sum-checked binary 
tape capable of bring immediately read into any program-specified storage location of the 
computer. The use of such a subroutine is obviously advantageous in a calculation where 
one has intermediate data which required a long time to compute and which must be reused 
several times. This routine is being submitted to the library of subroutines. 

A third program developed is a subroutine for generating spherical Bessel functions 
of imaginary argument, namely, 


‘ n * X * "-fix 'n + 1/2 *** * nd K n + 1/Z M where n = 0(1) 20. 


A fourth and closely related subroutine essentially generates the functions A n (x) 
and B r (x), n ■ 0(1)20, which are useful in two-center molecular integral calculations and 
are only sparsely tabulated. These are defined by 

oo 

A n (x) * J t n e” Xl dt 
l 

B„(x). /, n .-«*dt. 2(-l) n 5l_ [i 0 («)] 

-1 dx 

In addition, two more programs have been written and are currently being tested. 
The first routine consists of a complete rewrite of the previous matrix-element generator 
program. The second program, which is in production style for other users, is for the 
computation of single and double integrations (Simpson's rule) of numerically given func¬ 
tions and is a basic part of one procedure of evaluating multicenter integrals. 

All of these routines have been developed by F. J. Corbatd of the MIT Physics 
Department. 

179 TRANSIENT TEMPERATURE Of A BOX-TYPE BEAM 

The transient temperature and stress response of an extruded rectangular cross 
section exposed to a step function heat-flux input has been calculated by L. A. Schmit of 
the MIT Aero-Elastic and Structures Research Laboratory making use of an uncoupled 
finite-difference procedure. F.r r idiation and convection losses were neglected, and the 
thermal properties of the structure were assumed constant. The calculations have been 
compared with laboratory experimental results for the first 10 seconds of exposure, and 
the agreement is considered satisfactory. The experimental results referred to here were 
obtained by Lt. F. L. Williams, USAF, and are reported in his Master's thesis for the 
Department of Aeronautical Engineering. The title of Lt. Williams' thesis is "The Com¬ 
bined Effects of High Intensity Heating and Dynamic Loading on a One Cell Box Beam." 

It should be noted that there is a chordwise variation in the heat flux. This is due 
to laboratory furnace edge effects. This chordwise variation is taken into account in the 
finite-difference solution carried out on the Whi rl wind computer. An approximate analytical 
solution for the temperature response has been obtained for the special case where the 
heat flux is uniform with respect to chord and time. In the particular case under investi¬ 
gation the chordwise variation of heat flux is small, and the results obtained from the 
approximate analytical solution are in substantial agreement with the results of the finite- 
difference solution. 
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I HO CROSSCORRELATION OF BLAST-FURNACE INPUT-OUTPUT DATA 

Thin problem was described in Summary Report No. 38. The combined program 
(or computing and plotting the correlation function and its Fourier transform has been 
completed. An estimated 2 hours of computer time will be required to complete test runs, 
after which a summary of operating instructions will be compiled for the use of program¬ 
mers wishing to use the program. 

This work will successfully complete the project initiated under Problem 108. A 
further phase of the research on application of advanced statistical techniques to operations 
research problems is now in its preliminary stages, and a future problem may be initialed 
to carry out the computational phases of the work. 

This problem has been carried out by Richard G. Mills, formerly studying in the 
Department of Electrical Engineering at MIT. 

183 BLAST RESPONSE OF AIRCRAFT 

The results of numerical integration of the blast-response problem by Milne's 
method were found unsatisfactory because of the existence of a discontinuity in accelera¬ 
tion. The formulation was then revised using the Runge-Kutta method. The time history 
obtained agrees very well with the corresponding case computed by desk calculators ac¬ 
cording to a closed-form solution for the quasi-steady damping case. 

A parametric study of peak response was made, totaling approximately 200 cases. 
The basic tape was then modified to provide an automatic increase of peak gust ratio to 
determine the lethal gust for a given set of airplane and operation parameters. The com¬ 
puting by WWI for determining the lethal criterion of two typical airplanes is about to be 
completed. 

This study is being carried out for his doctorate thesis by H. Lin of the MIT Aero- 
Elastic and Structures Research Laboratory with the assistance of Y. Shulman. 

184 SCATTERING OF ELECTRONS FROM HYDROGEN 

The calculations described in the previous Summary Reports have been completed. 
These amounted to the evaluation of various two- and three-dimensional integrals which 
arose in the calculation of the differential amplitudes for elastic scattering of low-energy 
electrons from monatomic hydrogen taking into account polarization and exchange. The 
result* have been presented in the thesis of Maurice Newstein which has been submitted to 
the Physics Department at MIT in partial fulfillment of the requirements for the degree of 
Ph. D. They will also appear as a final report for the Committee on Machine Methods of 
Computation at MIT. Also it is planned to submit portions of the thesis to be published as 
an article in the Physical Review. 
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189 DISTRIBUTION OF GUSTINESS IN THE FREE ATMOSPHERE 

The first probability distribution of guiltiness can be obtained from the spectrum of 
the radar signal returned by clouds or precipitation. These meteorological developments 
are being used essentially as tracers of the atmospheric motions. 

The spectrum was computed from the autocorrelation function through a program 
devised by Douglas Ross. Six examples of such spectra were obtained, all of which were 
entirely satisfactory. 

This study is being carried out by A. Fleischer of the MIT Project on Weather 
Radar Research. 

190 ZEEMAN AND STARK EFFECT IN POSITRONIUM 

The computations in this problem, reported as being completed in Summary Report 
No. 39, actually extended over into the beginning of this period. A detailed report was 
given in Summary Report No. 39 by H. W. Kendall of the MIT Radioactivity Group. 

191 EARTHQUAKE EPICENTER LOCATION BY GEIGER'S METHOD 

Geiger's method of locating an earthquake uses arrival times of elastic waves at 
three or more observing stations. The method is basically a least-square process of fitting 
observed data to previously obtained empirical travel-time curves. This method has been 
mechanized, for an example using eight observing stations, as a part of D. R. Grine's 
S. M. thesis (Department of Geology and Geophysics, September 1954). 

The program iterates Geiger's method until corrections to earthquake location are 
smaller than inherent errors in the observed data. 

A generalization of the program is to be written for publication at a later date. 

194 AN AUGMENTED PLANE WAVE METHOD AS APPLIED TO SODIUM 

Summary Report No. 38 outlined the problem of reducing the order of the secular 
equation for energy levels in solids by first making proper linear combinations of the wave 
functions so that they are broken down into irreducible subspaces. 

Mr. M. M. Saffren of the Solid State and Molecular Theory Group at MIT haswritten 
and partially tested a program which chooses these linear combinations of wave functions 
as a function of the reduced wave vector of these wave functions and the Irreducible sub¬ 
space to which these linear combinations are to belong. When wave functions belong to the 
same order of neighbor in reciprocal space and are made into several linear combinations 
there is always a possibility that not all these combinations are linearly independent. This 
is extremely troublesome in the solution of the secular equation so the program has been 
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provided with the facility of disregarding linearly dependent solutions. The program has 
the additional feature of being able to choose the linear combination with the least number 
of nonzero coefficients if all that is desired is a linear combination which takes functions 
into an irreducible subspace without consideration of how basis functions are chosen in that 
space. This feature is a valuable one when the time for formation of the matrix elements 
between these linear combinations is strongly dependent on the number of nonzero members 
in the combinations. 

The program operates by picking up a wave vector from magnetic-core memory 
and then finding the subgroup which belongs to this wave vector. The subgroup is coded in 
three WWI words. The program then picks up the simple code words which tell into what 
basis function of which irreducible subspace the linear combination of wave functions is to 
go. If no basis function is specified the program chooses the minimal linear combination 
for the subspace. The result is six WWI words which are then fed into a matrix-generation 
and diagonalization routine (developed by Dr. Howarth under Problem 147). Finally the 
program picks up a series of reciprocal lattice vectors each of which characterizes a 
different order of neighbor and generates all the different reciprocal lattice vectors for 
each order of neighbor. These vectors are then processed so that only reciprocal lattice 
vectors which cannot be obtained from one another by operations of the original wave-vector 
group appear. This ensures against a secular equation with rank lower than order. These 
reciprocal lattice vectors are then fed into the matrix-generation and diagonalization routine 
and this routine is entered. 

It should be mentioned that Mr. Saffren's routine can, at the moment, only be 
applied to solids which arc face- or body-centered crystals. This is because Dr. Howarth's 
routine applied only to such crystals. 

More details regarding the theory of the calculation can be found in the Quarterly 
Progress Reports. Solid State and Molecular Theory Group, MIT. 

195 INTESTINAL MOTILITY 

A study o£ th* effect of radi.llon upon the gastrointestinal tract of rabbits is being 
carried out for Dr. J. T. Farrar of the Gastroenterological Section of the Evans Memorial 
Hospital by Miss D. Hamilton of the MIT Servomechanisms Laboratory. 

During this period several runs were made, but the main effort has been to program 
a routine for plotting results and to use techniques which reduce truncation errors in the 
Fourier transformation. Since difficulties arose with the automatic scope routine, special 
scope routines are being developed. 

The autocorrelations and the Fourier transforms cited above use programs devel¬ 
oped under Problem 107. The scope routine for plotting and the techniques for reducing 
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truncation errors were developed unde r Problem 171. These problems have been described 
by Douglas Ross in previous reports under their respective problem numbers. 

199 LAMINAR BOUNDARY LAYER OF A STEADY, COMPRESSIBLE 
FLOW IN THE ENTRANCE REGION OF A TUBE 

In connection with research on heat-transfer coefficients, recovery factors, and 
fraction coefficients for supersonic flow of .air in a tube, a theoretical investigation of the 
characteristics of the laminar boundary layer in the entrance region has been carried out 
by Dr. T. Y. Toon^ of the MIT Mechanical Engineering Department. Partial differential 
equations of continuity, momentum, and energy were developed for the boundary layer. 
These were then transformed into a series of ordinary differential equations to be solved 
by WWI for several entrance Mach numbers and thermal conditions at the tube wall.. 

First, solutions are to be obtained for the case of constant viscosity and thermal 
conductivity. Then, the effects of temperature dependence of these properties are to be 
studied. 

Solutions of the second set of differential equations for the case of constant vis¬ 
cosity and thermal conducts uy have been obtained for six cases of different entrance 
Mach numbers and thermal conditions at the tube wall. This represents another 15% of 
the entire job to be done by WWI. 

Work is being done to solve the third set of differential equations corresponding to 
two of the six cases mentioned above. 

200 A STUDY OF RECURRENT EVENTS 

The investigation of various digital methods for identifying one of two possible 
Bernouilli sequences has been completed. This work has been done by Miss B. Jensen for 
Dr. G. P. Dineen, both of the Lincoln Laboratory at MIT. 

Although various methods were tried in order to determine the most efficient meth¬ 
od, no rigorous optimization scheme was carried out. The use of "success runs" 1 and 
their recurrence times as a means of distinguishing a Bernouilli sequence with high prob¬ 
ability has been investigated using computer-generated sequences. A similar study for 
moving average detectors was also carried out. The use of Whirlwind I made possible the 
acquisition of a large amount of numerical data on these detectors. 


1. Feller. W. . Probability Theory and Its Applications, Wiley, New York. 1950. 
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201 STUDY OF THE AMMONIA MOLECULE 

Since the last report the self-consistent field solutions for the ammonia molecule 
have been computed. From these, orbitals which represent the ground state of the mole¬ 
cule have been constructed. These, in addition to the excited wave functions for the mol¬ 
ecule. will be used to get the best approxlmution to the wave function of this molecule within 
the approximations used. The interaction of the lowest state and the excited state will give 
rise to a matrix of interaction. The diugonalization of this matrix will be carried out on 
Whirlwind through the use of a prog ram written by Alvin Meckler and improved by F. Cor- 
bato of the Solid State and Molecular Theory Croup at MIT. The calculation of the matrix 
of interaction is now in progress. 

203 RESPONSE OF A FIVE-STORY FRAME BUILDING UNDER DYNAMIC LOADING 

The five-story building has been approximated by a system of five mass points 
(representing the floors) connected by nonlinear acting springs (simulating the building 
columns). This representation of a building frame is conventional for frames with stiff 
floors so that any failure will occur in the columns. For frames up to five stories high it 
seems to be a good approximation and results in equations which are readily adaptable to 
numerical calculation. 

The equation describing the motion of a mass point is: 

m .v 5 “„>= p <‘ n ) \ 

where 

m g = mass s 

y s = horizontal displacement of mass s 

y = acceleration of mass s 

P(t n ) = load at time t^ (impulse) 

Y s * loading constant for mass s 

R s + I = resistance of spring above mass s 

R s = resistance of spring below mass s 

The loading is calculated from pressures determined from the actual analysis of a 
five-story building under dynamic loading. 

To find the value of the resistance of any spring s at any time t it is necessary to 
calculate an assumed elastic value of the resistance and then compare this value with the 
maximum possible spring resistance at this time. These maximum values of the spring 
resistance at various times are known in advance and are included in the program as a 


series of constants. The resistance calculation is based on the previous value of the 
resistance and the change in the relative displacements of the two masses adjacent to the 
spring. Therefore: 

R 1 (t ) ; R (t .) t K ,1a ,(l ) - a ,(l ,)1 

» ' n' s'n-r sl[ s,s-l'n' s.s-l'n-rj 

Ir *(l ) I <R , .(l ) 

I * ' n' I s(max)' n' 

where 

R s *(t n ) s value of the resistance at time (t^) 

R s (t n j) = previous value of the resistance 

K s j = slope of the equation relating a and R in the linear range (for spring ») 

a s s l^ l n^ = re ^ al * vc Reflection of mass s with respect to mass s-1 at time (i^) 

R / v (t ) = maximum value of the spring resistance of spring s at time t 

s(max)' n ’ ° n 

This procedure allows for increasing or decreasing relative deflections which may result 

in either positive or negative values of the resistance. One disadvantage of the procedure 

is that errors in the value of R s (* n ) tend to accumulate. For this problem this disadvantage 

will not result in obvious errors. 

For the integration of the five nonlinear second-order differential equations, the 
second difference equation, 




To check the program, simplifications were introduced in the original program and 
the response of a three-story frame building was determined. These results checked very 
well with a hand calculation and also compared favorably with experimental results. The 
five-story frame has been analyzed under two similar loading conditions and has produced 
results which agree rather closely with a very elaborate hand calculation. 

The program is now being modified slightly to include the effect of vertical load 
plus dead load. It is hoped that this refinement will help in producing results even closer 
to those achieved in the hand analysis. The handling of the loading constants is also being 
modified to smooth out the data, particularly in the neighborhood of sharp changes in slope. 
The modified program will then be used to determine the response of the five-story building 
for about 30 different loading conditions. 

This study is being carried out by C. W. Johnson of the MIT Department of Civil 
and Sanitary Engineering. 
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204 EXCHANGE INTEGRALS BETWEEN REAL SLATER ORBITALS 

The purpose of Problem 204 is the computation of exchange integrals between Slater 
atomic orbitals located on two separate centers. The program designed to accomplish this 
can, with slight modifications, be generalised to give two-center hybrid and Coulomb inte¬ 
grals as well as exchange integrals. The efficiency with which these extra tasks can be 
accomplished by use of the program is dubious. However, this generalization has been 
incorporated. The program has been coded in its entirety. The tapes have been prepared 
and the testing procedure initiated. 

The computation of one of the auxiliary functions, 

I(n./3, x)» J yV* y dy. 

X 

is expected to introduce additional difficulties resulting from cancellation errors in the 
generation of the functions. Since P. Merryman has discovered that this function can be 
evaluated (in the desired range) as the quotient of two series of positive terms, this diffi- 
culty--il It docs appear--is not expected to prove insurmountable. 

This work is being carried out by P. Merryman under the direction of Professors 
Mulliken and Roothaan of the University of Chicago Physics Department. 

208 INTERCEPTOR FLIGHT-CONTROL PROBLEM 

This problem on interceptor flight control was described tn the last Summary Report 
as a system of 16 nonlinear differential equations and was solved using the Gill method, an 
adaptation of the Runge-Kutta fourth-order procedure. 

The system was varied in this quarter, a new equation substituted, and a constant 
changed to give a system of 1} differential equations. Two solutions of this system were 
run on Whirlwind using CS U. The stability of one solulion and the instability of the other 
demonstrated the sensitivity of the whole system. 

This represents the end of the digital computer phase of this problem. The solution 
was carried out by M. Merwin of the MIT Dynamic Analysis and Control Laboratory. 

209 NUMERICAL SOLUTION OF HOMOGENEOUS LINEAR DIFFERENTIAL 

EQUATIONS WITH QUADRATIC POLYNOMIAL COEFFICIENTS 

This problem was selected by Dr. J. C. P. Miller of Cambridge University, Cam¬ 
bridge, England, as one that not only would provide a generally useful and interesting 
routine but could be completed during his brief visit at the Digital Computer Laboratory. 

The method of solution makes use of a recurrence relation for terms in the Taylor 
expansion. A detailed discussion may be found in Summary Report No. 39. 


During the quarter a few more test runs were completed, and the problem has ter¬ 
minated. 

210 RESIDUE INDICES AND PRIMITIVE ROOTS 

This problem was also selected by Dr. J. C. P. Miller of Cambridge University, 
Cambridge, England, to be coded on WWI during his brief visit. The problem is a famil¬ 
iar one: given aprimep and a small integer a, determine the least e such that a C -l is a 
multiple of p. 

The numerical procedure selected was to determine from r. s a* mod p, r. + j = ar^. 
This was then divided by p to give This was tested for fl and the process 

repeated, counting the number of repetitions. 

Dr. Miller plans to code this problem for five or six computers as part of a com¬ 
parison test of speed and capabilities. The coding on WWI was completed before Dr. 
Miller's return to England. 

211 SERVO RESPONSE TO A COSINE PULSE 

This problem was concerned with extending the curves in Figure 25-16, p. 705 of 
Instrument Engineering , Volume II, by Draper, McKay, and Lees. Work on this problem 
has been completed. The results have been plotted and are being considered for inclusion 
in Volume III of the above-mentioned book. The original problem was to find numerically the 
maximum response x(t) for some parameter pairs (£, Tp/T^O.Ol £ £ 5 5, 0. 1 S T^/T^ 
£ 10 in the equation 

(l/w n *) >1 + (*£/w n ) x t x * G(t), 

where T » 2 if/w , it = dx/dt, V = d^x/dt^. and where G(t) e 1/2 [l-cos(2irt/T )] for 0 S t 
n n P 

S T , and G(t) « 0 for t > T^. As work progressed on this problem it was decided also to 

treat in the same manner the similar problem corresponding to 

G(t) • 1/2 [ I -cox (2rt/T p ) + 2t(T n /T p ) .in w p( ] 

for 0 s t 5 Tp, and G(t) = 0 for t >T^. 

Coding of the problem using the CS II system was accomplished by Dr. J. M. Stark 
of the MIT Instrumentation Laboratory. Application of the curves obtained is described in 
the above-mentioned reference. 
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212 DISPERSION CURVES FOR SEISMIC WAVES MULTILAYERED MEDIA 

A program hat* been written which calculated dispersion curves for surface seismic 
(Rayleigh) waves on multilayered media. The calculation follows that outlined by Haskell 
(••The Dispersion of Surface Waves on Multilayered Media, " N. A. Haskell, Bui. Seismo- 
logical Society of America, V. 43, No. 1, January 1955, pp 17-34). The only difference is 
that, rather than assuming a phase velocity and iterating to find a wave number, the reverse 
procedure is followed to avoid difficulties which would arise at a phase-velocity minimum. 
The program has been written for as many as five layers. When the main body of the pro¬ 
gram is working properly, a routine will be added to differentiate the results numerically 
in order to obtain the group velocities. 

This work is being carried out by K. Vozoff of the Department of Geology and Geo¬ 
physics in partial fulfillment of the requirements for the degree of Ph. D. 

213 INDUSTRIAL PROCESS CONTROL STUDIES 

The Machine Design Division of the Mechanical Engineering Department is engaged 
in determining the characteristics of a physical system or process from the random fluctu¬ 
ations in the input and output records during the normal course of operation of the system 
or process. The idea has been developed by T. P. Goodman, C. M. Chang, and J. B. 
Reswick from work done by N. Wiener, A. Tustin, and Y. W. Lee and has been found to be 
promising in its applications both to linear systems and to some nonlinear systems under 
actual operation in industry. It is believed to have ce rtain advantages over the conventional 
approaches (e.g., step and frequency responses). 

It may be shown that certain statistical functions of the random data (the auto- and 
crosscorrelations) bear the same relationship with each other as the random input and 
output themselves. Since these statistical functions possess many desirable properties 
that the random functions themselves do not have, they are used for the determination of 
the system characteristics. The computation of these statistical functions is hence a most 
essential step in the current investigation. Such computations involve a great number of 
successive multiplications and averaging and can be done most readily on an automatic 
computer. 

In the past 3 months, some exploratory work on the new technique was done. Sta¬ 
tistical functions on data taken from a real flow system and some artificial systems were 
computed. These functions, after analysis, gave system characteristics which compared 
very favorably with the system characteristics obtained by the conventional approach. The 
essential results are reported by C. M. Chang in his Mechanical Engineer's thesis, entitled 
"A New Technique of Determining System Characteristics from Normal Random Operating 
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Records." 

Efforts will be continued in much the same direction in the near future. Arrange¬ 
ments have been made with some industrial concerns so that data on actual industrial sys¬ 
tems can be taken in their plants. The statistical functions can then be computed and ana¬ 
lyzed and the system characteristics determined. It is hoped that with further experience 
and study many of the present difficulties can be solved and new light will be shed on the 
characteristics of industrial processes, 

214 INTERVAL DISTRIBUTION 

This problem arose in the textile industry. A suitable distribution for the "slubs” 
or thickenings on an artificial yarn was needed. The one suggested, which could be pro¬ 
duced by a simple mechanism, was this: Take a series of events occurring at regular in¬ 
tervals a and delay each event in a random manner, the delays being distributed about a 
certain mean in a gaussian way. 

The problem is to calculate the distribution of the intervals between successive 
events produced in this way. The method was the Monte Carlo one of making a long series 
of events in the prescribed manner and calculating the interval distribution. There was no 
loss in generality in taking the mean delay as aero (since it is irrelevant) and the standard 
deviation of the delays as 1. The only parameter is then a, the original interval. 

For large a, the interval distribution itself is normal, because the sequence of 
events is almost never disturbed by the random deviation. For small a, the distribution is 
like that of a Poisson distribution, except that for large intervals it decays as a gaussian 
distribution docs. The intervals are not uncorrelatcd, however. The interesting region to 
investigate is that around a= 1 in which the distribution changes from one kind to the other. 

To produce the random normal deviates, uniformly distributed random numbers 
were added. The difference equation a n ^j = a n + a n 1 ( moc * U produces such numbers. 
These series were calculated double length to reduce the dangers of a periodic solution. 

Twelve successive values were added, and, after subtracting the mean, which is 6, 
a reasonably good approximation to a normal deviate with standard deviation 1 should be ob¬ 
tained. This distribution was tested by distributing into cells of width 1/4 in the range -2 
to *2 four runs of 10,000 cases. A count of cases falling outside the range was kept as a 
check. The results were satisfactory. 

Having the spacing a and random deviates Vj v^ v^, the intervals, if the sequence 
were preserved, would be a ♦ v^ *Vj, a 4- v^ - v^, etc. However, some of these would ap¬ 
pear negative, because the sequence is altered. We therefore keep the last ten intervals in 
case the next event falls somewhere in this region. 

Suppose three successive intervals are ij, i^, ij and i^ is negative. Then trans- 
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posing the twoevents of which i^ is the Interval gives the new intervals i^ ♦ -i^. lj ♦ 

If the ten intervals kept have been "sorted" and are positive, when the next interval 
is found, if it is negative, the transposition process, carried down the series of Intervals, 
will eventually make them positive. If, however, the last interval stored were to turn out 
negative, the process would have failed because more than ten should have been stored. 

As a new interval is added to the series of ten, the oldest member of the series is 
added to the distribution and discarded. 

For this problem, cells of width l/4 between the ranges 0 and 8 were used, and, as 
a check, a count was kept of intervals greater than 8 and apparently negative intervals, 
though in fact none were found. 

Two runs for each of the values a = 1.0, a = 1.5, a = 2. 0 were done, with 10,000 
intervals in each. Before beginning to accumulate the distribution, 20 intervals were cal¬ 
culated to fill the store of intervals properly. 

The results were treated like experimental results and v hen plotted gave smooth 

curves. 

This study was carried out by D. W. Davies of the National Physical Laboratory, 
Teddington, England. 

215 PLANT SURVEYS BY STATISTICAL METHODS 

The application of automatic feedback control to existing industrial heat transfer 
processes requires knowledge of the transfer function of the process. Such transfer func¬ 
tions are ordinarily determined by observing the response of the process to a sinusoidal 
variation of its input variable. The aim of the present study by S. Margolis of MIT's Re¬ 
search Laboratory of Electronics is to perfect methods for determining transfer functions 
from formal operating records without introducing any variations in the input variable other 
than the random fluctuations already existing. 

A small heat exchanger installed in the Process Control Laboratory of the Electri¬ 
cal Engineering Department has been used for experimental studies. Random fluctuations 
of the input and output variables were recorded on a Sanborn Two-Channel Recorder. This 
data was reduced manually to numerical form. The autocorrelation function of the input and 
the crosscorrelation between input and output were computed by Whirlwind 1 for two separ¬ 
ate sets of data, using routines developed by D. T. Ross of the Servomechanisms Labora¬ 
tory. The transfer functions determined from the correlation functions agree with those 
determined by subjecting the system to sinusoidal variations. 

216 ULTRASONIC DELAY LINES 

This problem was undertaken in cooperation with the Arenberg Ultrasonic Labors- 
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lory because of the importance of the result, for MIT's Lincoln Laboratory. The problem 
is concerned with the optimum design of an ultrasonic delay line. 

For the optimum packing of microseconds delay per cubic inch, radial or "a.terl- 
ated" symmetry paths arc used. This design was developed by the workers of the Bureau 
of Aeronautics of the U.S. Navy. These designs have given way to MS or "Multiple Sym¬ 
metry" designs where many of the regular facets are tilted. 

A routine has been developed for W WI that will select by iteration an optimum set of 
radii for such designs. 

217 VARIATION-PERTURBATION OF ATOMIC WAVE FUNCTIONS AND ENERGIES 

The wave functions and energies for low-lying stationary states of helium through 
neon and their isoelectric series have been calculated approximately from a variational 
method. * The method consists of minimizing the integral 

E *fji H t#r dt (1) 

where H is the Hamiltonian of the system and t includes all spatial and spin coordinates. 
\j/, the trial wave function, is a properly normalized linear combination of determinantal 
functions whose components arc hydrogen-like. The nuclear charge in the hydrogen-like 
functions is replaced by variational parameters. 

The integration in (1) may be carried out analytically leaving E as an algebraicfunc- 
tion of the variational parameters and reducing the problem to one of finding the minimum 
of a rather complex expression. 

A program for Whirlwind I, which computes the best wave functions for any atom 
containing Is, 2p, and 2s orbitals, has almost been completed. It will probably be extended 
shortly to include the 3s and sp orbitals. The extremum for E is found by expanding it in 
a Taylor's series to second order about a point assumed to lie near the extreme point, dif¬ 
ferentiating the series with respect to the parameter increments, and solving the resulting 
series of linear inhomogeneous equations for the coordinates of the extreme point relative 
to the zero point of the Taylor expansion. 

After the initial variational treatment, the wave functions obtained can be improved 

by adding to them a linear combination of correction terms similar in structure but corre- 

2 

sponding to higher excited states. 

\lr «\lt * y c a 

T corrected T variational A c i *i 

1. Morse, P.M., Young, L.A., Haurwitz, E.S., Physical Review 48, 948 (1935) 

2. Yilmaz, H., Ph. D. Thesis (unpublished), MIT Physics Department (1954) 
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The c^'s arc determined by first-order perturbation theory. In cases where the 0^ contain 
parameters not already fixed by the first variation, a second variation involving the unfixed 
parameters is carried out. 

Several of these perturbation cor rcctions involving the (Is^)'s and (ls£p)^p, *p con¬ 
figurations of helium have been car ried out on Whirlwind by Arnold Tubisof the MIT Physics 
Department. The results are discussed in Quarterly Progress Report No. 14 of the Pro¬ 
ject on Machine Methods of Computation and Numerical Analysis, MIT, December 15, 1*)54. 
The work on these pc rturbation corrections is being continued at present. It has been found 
that both the variation and perturbation calculations may be broken down into well defined 
subroutines which greatly facilitate the programming. 

218 TRANSFORMATION OF INTEGRALS FOR DIATOMIC MOLECULES 
Eigenvalue equations are frequently encountered of the form 

[h] X * X [ s] X. ( 1 ) 

where [h| is He rmitian and £sj is real, symmetric, and positive definite. Standa rd meth¬ 
ods based upon successive orthogonal or unitary transformations can be applied to solve 
equations 

[H'jYsAlY. U) 

where I is the unit matrix. To convert equations (l) into this form by the Schmidt process, 
a lower triangular matrix ^Mj is constructed which has the property 

[m] [s] [m] . 1. (3) 

Here a tilde (~) denotes a transposed matrix. Applying transformation (3), equations (1) 
are reduced to the form (2) with 

M - MHffl <«) 

and X = [ 5 ] Y. (5) 

Whirlwind programs have been coded which construct [m J from a given positive- 
definite matrix [ s] of arbitrary dimension n x n and carry out the transformations (4) and 
(5) for any real symmetric n x n matrix ^hJ and n-vcctor Y. 

Given [sj . the Schmidt matrix [m] is constructed as a product of triangular ma- 


M *[>"]•■• K] 


The matrix 
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[s m ]. [b- 1 ]... [b 1 ] [s] [b 1 ] ...[B m -*] 


has the unit matrix I m * as its upper left-hand suLmatrix of orde- (m - 1)^. The remain¬ 
ing elements arc in general nonzero. The well-known Schmidt orthonnrmali/.atinn process, 
applied to the basis vectors of [s** 1 ] , normalizes the vector of index m and makes it orth¬ 
ogonal to all vectors with indices 1, .... m - 1. Vectors of index lower than m are not af¬ 
fected by[ b* 11 j . which is defined by 




i i, t < m or i, j > m 


B = -S B ; j < m 

mj mj mm 

m- 1 

B m = fs m - I (S m ) 2 l* l/ * 

mm L mrn mj J 

jo 

The independent elements of the symmetric matrix S are stored in lower triangle 
form, i. e. , in the order *11 S il S li *3 1 s u s 13' etc. The program which constructs ^ Mj 
is arranged to replace f sj by [Mj as the latter is computed, and only four additional num¬ 
ber registers arc needed. The programs for transformations (4) and (5) have similar prop¬ 
erties. 

This study is being carried out by R. K. Nesbct of the MIT Solid State andMolecular 
Theory Group. 


220 PROBLEM ARISING FROM AN ALGEBRA 

As a sample problem to help in his study of WWl, D. W. Davies of the National 
Physical Laboratory, Teddington, England, selected the study of the algebra whose basic 
elements consist of two associated patterns which can be combined as follows. 

Consider the pattern (a) representing one way of joining in pairs all of six spots. 
For six spots there are 15 such patterns. Taking two patterns, (a) and (b), if they are 
supe rimposed on one set of spots a number of loops are formed. This number will be called 
J?(a.b). and 1 < 9 S 3. 

In general the re can be <JN spots and (2N- 1)(2N- 3) . . . (3)( 1) ways of joining them in¬ 
to patterns. For any two patterns a, b, 1 S $ (a, b) < N. 

Suppose we take two fixed patt-rn.: a, band a variable pattern x. The numbers J?(a,x) 
and $ (b, x) can be used as coordinates to establish the pattern x in a cell of an N x N ma- 
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trix. Writing in each cell the number of element* with these coordinates, we have a distri¬ 
bution of all possible elements according to their values f (a, x) J (b, x). 

It turns out that the numbers in the matrix depend only on f{ (a, b) and not on the choice 
of a and b. For this reason, the matrices are symmetric, and there are N of them for the 
patte rns of degree N, corre sponding to Q (a. b) 1,2, ...N. 

The work involved in computing these matrices by hand becomes large after N = 3. 
Cases N 4. 5 are done quickly on Whirlwind, but the six matrices for N = 6 took 45 min¬ 
utes to compute. Case N = 7 was therefore not attempted. 

The computed matrices showed a pattern of zeros, which could be expressed by the 
two inequalities 

Jf U,») (b.x) S N f J? (a.b) 

f(a.x) - f(b,x) < N -jfU.b) 

A proof of these relations (which arc really identical) was obtained. 

By examining the matrices carefully, a relation between them was discovered. This 
relation enabled the matrix of order N for $ (a.b) n to generate the matrix of order N + 1 
for^(a,b) = n f 1. A proof of this relation was obtained. 

With this relation, all the matrices of order N t 1 can be obtained except the one 
JP(a.b) = 1. The work of calculating the matrices would have been greatly reduced if it had 
been known. 

The symmetry of the malt ice* $ (a. b) 1 and their invariance to the choice of a and 
b can easily be shown. This, however, docs not prove the invariance for the other matrices, 
and no proof has been found for it. 

221 COURSE 6. 25, 1954 

Fifty-nine students enrolled in course 6.25, entitled Machine-Aided Analysis, which 
was given in three sections during first semester, 1954, by Professors W. K. Linvill, 
R. C. Booton, and C. W. Adams. M. chines which aid in the solution of mathematical equa¬ 
tions encountered in engineering were studied. TAC, the three-address summer session 
computer, was used by the students as an exercise in digital-computer programming. A 
detailed description of TAC may be found in Summary Report No. 39. 

The following problem was solved: 

"Consider the Van dcr Pol Equation: 

V (y* - U y ♦ y = o, where 

V(0) = a 
y(0) = b 


For this equation assumed values of < , a, and b are chosen so that y 

'max 

< 10. y < 10, andy < 10. Solve the rquaiionusing intervals, h * 0. 1. 
max max 

Print out the values of y for integer values of t and carry out the solution 

until the value of y has passed its second positive peak. Scale factor the 

solution so that all of TAC's 8 digits are used. Use the rectangular rule of 

integration. " 

Each student prepared his own Flexowriter tape and wa* present at the running of his prob¬ 
lem. This problem took approximately i minutes of W WI computer time to solve. 

222 HELICOPTER ROTOR STABILITY 

The equations of motion of a flexible helicopter blade with constant geometric and 
mass properties in forward flight were developed. This analysis considered two degrees 
of freedom of the blade, one, rigid flapping, llie other, elastic bending. Previous analy¬ 
ses made of stiff blades have shown that the transient solution decreases in stability with 
increasing forward speed to tip speed ratio. The forced solution has always been assumed 
to be a Fourier scries expansion in the azimuth position of the blade. 

A preliminary program, conducted as a Master's thesis, has shown that inclusion 
of the bending degree of freedom results in a prediction of decrease in transient stability 
which is in much closer agreement with experimental results than was obtained earlier. 

The proposed program consisted of two parts 

1. To supplement the results alr»*!dy obtained in the preliminary program by desk 
calculators; 

2. To integrate the equations of motion directly, thus to gain an insight into the be¬ 
havior of the forced as well as thfc transient solution. The equations of motion are given by 
a pair of second-order ordinary differential equations. 

In part 1 the problem was to expand an 8 x 8 determinant whose elements arc poly¬ 
nomials in a, the exponential factor of the transient solution. The expansion results in a 
16th degree polynomial in o. whose most positive root determines the degree of stability 
of the solution. This was done for various sets of the parameters of the problem. 

The numerical method makes use of the polynomial-transformation theorem and the 
Routh-Hurwitz criterion for positive real roots. (Ref. Draper, McKay, and Lees - Instru- 
ment Engineering, McGraw-Hill, 1953, Chap. 21.) 

The results of the first part checked with the hand results, and additional results 
served to throw more light on the behavior, stability-wise, of rotors with high moments of 
inertia at high speeds. 

Part 2 of this program consisted of a numerical integration, using the fourth-order 
Kunge-Kuita method, of the two simultaneous, ordinary linear differential equations (*e- 
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cond-order) with periodic coefficients. 

The second part, whose purpose was to provide visual results of the flapping motion 
of helicopter rotors with different inertia properties at various speeds, will be completed 
during the next quarter. Twenty-five cases will be run. 

Conclusions that may be drawn from this program are: (1) Higher harmonics should 
be included in the assumed solution; and (2) Whirlwind can be used to advantage in solution 
of problems of this nature. 

This study is being carried out by Y. Shulman of the MIT Aeronautical Engineering 
Department. 

225 NEUTRON-DEUTERON SCATTERING 

The problem consists of the calculation of S wave phase shifts for neutron-deuteron 
scattering by means of a variational principle. Stationary expressions are obtained for the 
quartet (J = 3 /l) and doublet (J = 1/2) phase shifts, which fully determine the low-energy 
scattering. These expressions involve only known functions and a so-called "inside wave 
function" y, which is the difference between the actual wave function and its asymptotic be¬ 
havior. Any trial function for y which goes to zero when cither neutron is far away will 
automatically satisfy the boundary conditions of the problem. The advantage of this form- 
.ulation is that it enables us to take account of polarization, i c the deformation of the 
deutcron while all three particles are interacting. 

If we make an expansion of our stationary expression in powers on the energy, the 
first two terms will give the scattering lengths and effective ranges, in complete analogy to 
the two-body scattering. Results of interest will be the dependence of the scattering lengths 
and ranges on the exchange nature of the forces and their comparison with experimental 
values. 

The Gaussian well is used for the two-body potentials, and tensor forces are not in¬ 
cluded. The trial function used is essentially Gaussian in the interparticle distances, with 
three variational parameters, \ rj . The numerical work consists of calculating a large 
number of algebraic expressions and a few numerical integrations as functions of \ , ^ , and 
tj and looking for a stationary point of the variational expressions. Unfortunately, one can¬ 
not tell beforehand whether this will take the form of a maximum, minimum, or inflection 
point. 

A program to calculate the quartet scattering using CS II has been written, and the 
checking of its various parts is almost complete. The calculation will begin with \ and p 
fixed at their "nonpolarization" values, and the variational expression will be plotted as a 
function of rj . X and p will then be varied in the neighborhood of the optimum value of 1) . 
It appears that a single curve will take between 5 and 10 minutes of machine time. 
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This study is begin carried out by L. Sartori of the MIT Physics Department. 

228 evaluation of difference diffusion equation 

The subject involves the one-dimension diffusion equation in the simplest form 

lit - n> 

i x i <m 


applied to a region 0 < x < 1 with the following conditions: 

Initial condition: * 0 

Boundary condition: jl£ „ 0 at x • l 
d x 


1. ^ c 1 and sin wt 

2 . At = -1 

O X 

3. if/ - Q At = 1 

d x 


at x 


0 


Explicit and implicit recurrence formulas of the coarse-grid difference approximation are 
used to derive the transient solution to this equation at the above conditions. 

The general recurrence formula can be represented as 

k+ i + ro k.i - k.ij-i, k +J > * *j. k ♦ - 9 ><*j + i. k - 2 *j. k + k>- 

By comparing the exact solutions o£ the continuous equation and the difference approxima¬ 
tion. a method is derived to optimise the conformation of the approximate transient solution 
to the true solution. The problem is to evaluate the merits of such a method over the 
conventional method. A rational way to compare these methods is to carry out complete 
solutions by digital computer for a period of time and to compare them with the continuous 
solution over the same period of time. 

With the exception of the interpolation process, the explicit ease of boundary con¬ 
ditions 1 and 2, above, are completed. The explicit case ol boundary condition 3 has been 
started. It is expected that the computation for the explicit case will be completed during 
the next quarter. 

This project is being carried out by A. T. Ling under the supervision of Professor 
S. H. Crandall of the MIT Mechanical Engineering Department. 


229 ROTATING CONTACT SEAL 

This problem seeks to evaluate the rubbing contact temperature in various designs 
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of a rotating contact seal for a gas turbine. The heat input at any part of the boundary in 
question is expressed as a linear function of the temperature at all other parts of the bound¬ 
ary, leading to a set of simultaneous algebraic equations. 

Two routines were written, and each ran successfully. Each routine generated 18 
matrices of order 15 x 15 and solved cor responding sets of simultaneous equations byCrout's 
method. 

This problem is being carried out by J. M. Bonneville of the MIT Mechanical En¬ 
gineering Department, and the results will be included in his Sc. D. thesis. 

230 BRIDGE ANALYSIS 

Saul Namyet of the MIT Civil and Sanitary Engineering Department is analyzing the 
dynamic response of a large number of bridges which vary in span and type. The bridges 
are represented by a one-mass system. This system is acted upon by a combination of 
time-variant external forces and a resistance function which is considered to vary with dis¬ 
placement in a predetermined manner. It is required to obtain the magnitude of a parameter 
(which determines the magnitude and time duration of the external forces) for which a pre¬ 
scribed displacement is achieved. The complete problem will consist of the analysis of the 
whole set of bridges for four different resistance functions. 

The solution of the problem consists of the solution of the following differential equa¬ 
tion: 

mil (t ) s P - R , 

' n' n n 

x(t n ) = the acceleration at time (t^) 

P s load at time (t ) 
n ' n' 

R n * spring resistance at time (t^). 

The equation is solved by use of the second difference equation: 

(At)- 2 >£ (t n )= *(. n+1 ) - z*(t n ) +*<!„_,). 

The program has been checked by comparing a limited number of response calcula- 
tions with hand solutions and graphical computations. Since the results have shown satis¬ 
factory agreement, the response is now being calculated for the complete set of bridges 
utilizing the first resistance function. 
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231 REACTOR RUNAWAY PREVENTION 

In this problem, M. Troost of the MIT Chemical Engineering Department has been 
attempting to detc rmlnc the safest reactor type and si/.e for MIT's proposed nuclea r reactor . 

The analysis involves a step-by-step algebraic calculation wjth an ite rative correct¬ 
ing routine. Many interesting results have already been obtained. 

232 ENERGY LEVELS IN A SPHEROIDAL POTENTIAL 

An investigation of the structure of heavy nuclei based on the Bohr-Mottelson model * 
is being carried out 4 by Kurt Gottfried of the MIT Laboratory for Nuclear Science. 

The first step in this program is the determination of energy levels and eigenfunctions 
associated with the motion of particles in a force field whose equipotentials are clHpuoids 
of revolution. Most of the required matrices, whose ranks range from four to fourteen, 
have now been diagonalized using F. J. Corbato's routine developed under Problem 172. 


1. Dan. Mat. Fys. Medd. 27, No. 16 (1953). 

2. Gottfried, K. and Weisskopf, V. F. , Quarterly Progress Report Laboratory for Nuclear 
Science, November 1954. 
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3. OPERATION 


OF WHIRLWIND I 


3.1 Summary 

Greater emphasis is now being placed on reliability and maintenance in order to 
provide the best possible performance. 

Reliability of the WWI computer system was slightly higher this quarter than in 
previous quarters. Applications time as reported by computer operators was 95% usable 
as compared with an average of 9 3% since September 1953. An account of the assignment 
of computer time to the various major activities is given in Fig. 3-1. 




Activity 

October [ November ] Der ember 

1-7 

• - 14 11-21 

22-24 

2*-4 

4-11 

12-14 

IV.24 

26-2 

t-V 

10-16 

17-21 

24->0 

Total 

Hour* 

Marginal Clinking 

4 

« 

1 

1 

1 

b 

7 

4 

b 

7 

6 

a 

4 

79 

Installation 

T 

0 


7 

7 

0 

7 

0 

7 

7 

7 

0 

0 

4V 

MaintvnaiK • 

21 


IH 

U 

If 

i 1 

14 

III 

24 

et 

10 

u 

10 

214 

Terminal Equipment 

Tealkng 

»4 

>6 

14 

21 

22 

21 

17 

11 

20 

n 

14 

u 

tl 

242 

Technician Instruction 

4 

4 

4 

* 

4 

0 

1 

" 

4 

i 

4 

0 

4 

40 

Scientific and Engineering 
Computitlon 

41 

17 

41 

II 

11 

11 

24 

2R 

to 

10 

24 

42 

21 

426 

Other Application* 

27 

)4 

19 

1b 

14 

64 

tj 

74 

9b 

16 

64 

44 

41 

714 

Total Hour* 

I4f> 

124 

1*1 

110 

14V 

111 | HO 

K. 

147 

144 


140 

41 

1824 


Fig. 3-1. Allocation of Computer Time 

3.2 Computer and Terminal Equipment 

During the past 6 months studies have been undertaken not only to improve system 
performance but also to train new personnel and to compensate for the recent transfer of 
experienced engineers to other areas. Specific phases of the present program are des¬ 
cribed below: 

1. A biweekly analysis of all equipment failures was initiated to point out the types 
of interruptions that occur most frequently and those that result in greatest loss of time. 


Periodic discussions of these records arid of unique »y stem-cngincc ring problems have 
helped to develop the single-system concept and make all personnel aware of the relation¬ 
ships of the various independent efforts to system performance as a whole. 

Z. As in the past, design weaknesses that show up are being corrected where 
practicable. Such weaknesses have already been eliminated in the core-memory sense 
amplifiers, in some of the buffer-drum circuits, in the marginal-checking control, and in 
the equipment-cooling system. 

i. New computer test programs have been written to facilitate routine marginal 
checking of some of the newly installed terminal equipment. However, power distribution 
for marginal checking must still be evaluated along with the test programs, and both may 
require changes in order to obtain satisfactory checking procedures. 

4. An extended training course in computer-system details is being conducted for 
new engineers and technicians for about J hours each week. Experience has shown th.it a 
thorough understanding of system logic is essential in order to locate trouble effectively. 
The course is supplemented by actual practice in locating simulated compute r failures. 

Over-all system reliability improved according to the daily estimates of computer 
operators. Their records show that 95% of the assigned applications time was usable, 
compared with the 93% average since September 1953. 

The biweekly interrupting-failure analyses initiated about Z5 September provide 
preliminary averages to supplement the figures above: 

1. Since Z 5 September, the average time between interrupting failures has been 
about 9.5 hours. (All system failures except those attributable to new instaIlations or 
modifications are included whether or not applications time was lost.) 

Z. During the 3-month period, the single category of failures causing most frequent 
stoppages was transient computer alarms which could not be immediately associated with 
faulty components or other system damage or deterioration. These alarms, over half of 
which arc not yet explained, accounted for 34% of the interrupting failures. 

3. During each biweekly period a few isolated failures, generally nonrepeating, 
caused major loss of time. Over the past 3 months, 10% of the interrupting failures have 
accounted for about 50% of the total time lost. 

It seems reasonable to expect some reduction in lost time as maintenance proce¬ 
dures are improved andoperating margins are widenedin the newer sections of the system. 
However, a large increase in the average time between failures seems doubtful, since 
detection and elimination of causes of transient errors are difficult. 
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4. CIRCUITS AND COMPONENTS 


4.1 Vacuum Tubes 

•4.1.1 Vacuum-Tube Life 

During the fourth quarter of 1954 the WWI computer operated for 1760 hours. Dur¬ 
ing 1954 total hours of operation were 6900, about 79% of the maximum possible hours in 
the year based on 365 twenty-four-hour days. This over-all utility factor is almost identi¬ 
cal to that for the year 1953. 

The tube complement of WWI, including the computer proper, magnetic-core me m - 
ory, magnetic drums, and miscellaneous input-output equipment, is now approximately 
12, 500. 

Vacuum-tube life has been calculated for the six most numerous types in WWI. 
These six types comprise about 80% of the total tube complement. A summary of this in¬ 
formation together with previous data is shown below. 

FAILURE RATE, PER CENT PER 1000 HOURS 


Tube Type 

First 

Quarter 

1954 

Second 

Quarter 

1954 

Third 

Quarter 

1954 

Fourth 

Quarter 

1954 

Quantity 
Now in 
Service 

7AD7/SR 1407/6145 

1.75 

1.6 

1.4 

1.3 

4237 

7AK7 

0.5 

0.6 

0.3 

0. 35 

3340 

6080/6080 WA/ 6AS7G 

1.1 

1.5 

1.4 

0.40 

718 

5965 

0.4 

0.3 

0 . 3 

0 . 26 

878 

6BL7GT 

0 . 3 

0.5 

0.5 

0 35 

494 


ANNUAL AVER AGE F AlLUkE RATE, PER C ENT PER 1000 HOURS 


Tube Type 

1952 

1953 

1954 

7AD7/SR1407/6145 

2.0 

3.3 

1.5 

7AK7 

0 . 26 

0.43 

0.44 

6080/6080 WA/6AS7G 



1 . 1 

5965 



0. 32 

6BL7GT 



0.41 


The failure rates for the 7AD7/SR1407/6145 and 7AK7 continue to be essentially 
constant. The failure rates for the 6080/6080WA. 5965, and 6BL7GT, which have been in 


service in quantity only about 15 months, continue to be encouragingly low, especially the 
latter two. 

In Figure 4-1 are given the vacuum-tube failures for the past quarter. A word of 
explanation is in order on several of the tube types The excessively high number of early- 
life failures of the 2D21 thyratron is dae to improper circuit design In .* relay driving 
application, the control and shield grid voltages were excessively negative during conduc¬ 
tion. This resulted in a high ion current to these electrodes, causing sputtering and rapid 
tube deterioration. Subsequent to a change in circuitry to reduce the magnitude of these 
voltages, no difficulty has been experienced. 

The 5963 is now the sixth most numerous tube in WWI, its quantity being only slightly 
less than that of the 6BL7GT. The failure rate for 5963's continues to be quite high, on 
the order of 4% per 1000 hours. Practically all of the failures have been electrical and of 
the nature of a change in characteristics, notably a decrease in plate current at zero bias. 
A detailed investigation has shown cathode interface impedance to be the cause of this low 
plate current in a considerable percentage of the failures. Previous life tests conducted 
here had indicated that this problem would develop (see Summary Report No. 33). Dis¬ 
cussions with the manufacturer indicate that steps have been taken to eliminate the problem. 

In Figure 4-2 are shown vacuum-tube failures for the year 1954. 

Failure rates for the years 1948 through 1954 for the 7AD7/SR 1407/6 145 and 7AK7 
arc shown in Figure 4-3. The points plotted arc for eachquarter of the years 1952 through 
1954 and the average of the years 1949 through 1951. The relatively high failure rates for 
both types during the first quarter of 1953 are a result of an intensive maintenance program 
undergone by WWI at that time. The sharp decrease in the 7AD7 failure rate in the first 
quarter of 1952 was due to a qhange in the marginal checking of flip-flops. The failure 
rate for the 7AK7 appears to he fairly constant at about 0.4% per 1000 hours. That of the 
7AD7/SR 1407/6145 appears to be tapering off between 1.0 and 1.5% per 1000 hours. 

In Figure 4-4 is shown ope ration time at failure for the 7AK7 and 7AD7/SR1407/6145 
for the years 1949 through 1954 normalized to the number of failures during the first 1000 
hours of operation. In both cases the number of tubes failing during a given 1000-hour 
period after about 10,000 hours of operation is only about 10-20% of the quantity failing 
during the first 1000 hours of operation. That the curve of the 7AD7/SR1407/6145 is smooth¬ 
er than that of the 7AK7 may be due to the fact that the total failures during the first 1000 
hours of operation are 249 and 33, respectively, and a small statistical v* riation represents 
a much larger percentage of the first 1000-hour quantity in the case of the 7AK7. 

4.1.2 Vacuum-Tube Research 

The life test of tungsten-nickel A31 Cathaloy Z-2177's has reached 2800 hours with 
one section conducting and the other cut off. There are no major changes from the last 
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Fig. 4-i. Failure Rate in WW1 
Tube Types 7AD7, SR-1407 61 IS and 7AK7 
19-18-1954 
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reading, considerable interface impedance remaining present on the conducting sections 
and much less on the cutoff sections. 

The life test of tungsten-nicke 1 All Cathaloy 5687‘s has reached 4200 hours with 
one section conducting and the other cut off Grid emission on the cutoff sections continues 
to increase; the grid currents on the conducting sections remain low. Only the same three 
of the ten conducting sections previously noted (Summary Report No. 19) show any cathode 
interface impedance, there is none on the cutoff sections 

As a result of interface-impedance difficulties encountered on the 5961 as mentioned 
in Section 4.1.1, a life test of new tubes was instituted during the past quarter. At 800 
hours the data is inconclusive. 

The 5965 life test has reached 6200 hours. Plate currents are holding up well; 
there is, however, considerable grid emission on the cutoff sections. 

A vacuum-tube pulse-chaructenstic tester is under construction for the tube lab. 
It will replace the existing breadboarded model which is much less versatile. 

4.2 Component Replacements 

Figure 4-5 lists the replacements of components other than tubes during the fourth 
qua rter of 1954 . 
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Fig. 4-4. Hours of Operation at Failure 
Tube Types 7AD7/SR-1407/6145 and 7AK7 
1949-1954 
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Total in No. of 

Service Failures 


0.0l- K f mica I0T. 5622 


Hours of 
Ope ration 


2000 - 3000 


0. 01- H f fixed 
ceramic disc 


1NI4A 



■I low R , 5 shorted 
Low R 

i> 


* b 
High R f 

l shorted, 1 low K, 


3AG (Buss) 


6100-ohm 
8-watt ♦ 5% 


9000-ohm 
1 / ?-watt ♦ 1% 


0. 45-megohm 950 

l/2-w.itt ♦ 1% 


0. 125-megohm 950 

1/2-watt ♦ 1% 


2500-ohm 
1/2-watt 
♦ 1% (Nobleloy) 


Sensitive 
(Mic roswitch) 

5:1 pulse, 
hi per all core 

7-pin miniature 



Out of tolerance 
Out of tolerance 
Out of tolerance 
Out of tolerance 

Out of tolerance 
Out of tolerance 
Out of tolerance 
Out of tolerance 

Out of tolerance 
Out of tolerance 
Out of tolerance 
Out of tolerance 

Out of tolerance 
Out of tolerance 
Out of tolerance 


Switch jammed 


Shorted turns in 
prima ry 


Fig. 4-5. WWI Component Failures October 1 - December 31. 1954 


5. ACADEMIC PROGRAM 


5.1 Programming Course 

The DCL CS 11 programming course was given once during this quarter. The course 
includes the following topics: relative addresses, temporary storage, floating addresses, 
preset parameters, programmed arithmetic, cycle counters, buffer storage, automatic 
output, post mortrms, and multipass conversion. The trxt for the course is a program¬ 
mer's manual written by staff members of the SfcEC Group. 

The 32 students enrolled during this quarter represented the following groups: 
Aero-Elastic and Structures Research Laboratory. Instrumentation Laboratory, Geophys¬ 
ical Analysis Group, Servomechanisms Laboratory, Meteorology Department, Physics 
Department. Electrical Engineering Department. Chemical Engineering Department. Civil 
Engineering Department, Mathematics Department, the School of Industrial Management, 
and Lincoln Laboratory. 


5.2 MIT Course 

The principal course on machine computation being offered at MIT in the fall of 
1954 was 6.25, Machine-Aided Analysis, a survey of computing techniques aimed largely 
at seniors in Electrical Engineering. This subject, first offered in the spring of 1954, had 
a fall-term enrollment of 59 students. It consequently was divided into three 20-man sec¬ 
tions taught separately by Professors Linvill, Booton, and Adams. Practice problems 
were planned to allow each student to use both a REAC and the Whirlwind I computer (sim¬ 
ulating the hypothetical TAC which had been developed for the 1954 summer session). The 
description of the problem done on TAC will be found in the report on Problem 221. 


5.3 Visitors 

Tours of the WWI installation include a showing of the film "Making Electrons 
Count," a computer demonstration, and an informal discussion of the major computer 
components. During the past quarter eight groups visited the computer installation. In¬ 
cluded in these groups were. 
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October 1 
November 3 

November 30 
December 14 

December 15 


Industrial Liiiion Office Personnel 

Professor Crandall's class (MIT Mechanical Engineering Depart¬ 
ment 

Emmanuel College Seniors 

Professor D. P. Campbell s class (MIT Electrical Engineering 
Department) 

A1EE-1RE MIT Student Branch 


6. APPENDIX 


6.1 Publications 

Project Whirlwind technical reports and memoranda arc routinely distributed to 
only a restricted group known to have a particular interest in the Project, and to ASTLA 
(Armed Services Technical Information Agency) Document Service Center. Knott Building, 
Dayton, Ohio. Requests for copies of individual reports should be made to ASTLA. 


The 

following is a list of memoranda 

published by the Scientific 

and Engineering 

Computations Group during the past quarter. 


• 

No. 

Title 

Date 

Author 

SR- 39 

Summary Report No. 39, Third 

Quarter 1954 


DCL-17 

A Program for Checking the 
Contents of Test Storage 

10-4-54 

A. Siegel 

DCL-22 

Utility Control Program 

11-22-54 

F. C. Helwig 

DCL-23 

Memo on SfcEC Performance 
Request 

11-19-54 

F. C. Helwig 

DCL-2* 

Director Tapes 

11-24-54 

F. C. Helwig 

DCL-25 

CS II Instruction STOP 

11-24-54 

F. C. Helwig 


58 


59 











